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ABSTRACT 

The dynamics of planetesimals in the circumbinary debris disc of the quadruple star system 
HD 98800 are investigated. Evolving a spherical shell of test particles from a million years 
ago to the present day indicates that both coplanar and retrograde warped discs could exist, 
as well as a high inclination halo of material. Significant gaps are seen in the discs, as well 
as unexpected regions of stability due to the retrograde nature of the stellar orbits. Despite a 
viewing angle almost perpendicular to the direction of the warp of the planetesimal disc it is 
still intersected by the line of sight for eccentricities of the outer orbit of 0.5 or less. 

Key words: celestial mechanics - planetary systems - methods: A^-body simulations - stars: 
individual: HD 98800 



1 INTRODUCTION 

HD 98800 is an interesting and unusual system of four 10 Myr 
old post T-Tauri K stars - two spectroscopic binaries A and B 
in orb it about one another - located in the TW Hydrae associ- 
ation jTorres et al] 1 19951 iKastner et al] Il997h . It has a large in- 
frared excess attributed to a circumbinary disc around the B pair 
jKoemer et al.l2000l : |Pra"to et alj200ll:lFurlan et alj2007h . Substan- 
tial extinction towards this pair is evidence that they are observed 
through some of this material jSoderblom et alJI 19981 pTokovininl 
1 19991) . A photometric variabili ty is also seen, but with no definite 
period jSoderblom et alj|l998|). Both the a pparent absence of CO 
molecular gas in the disc jPent et alj|2005i) and infrared spectrum 
modelling indicate that this is a T-Tauri transition disc that is just 
reaching the debris disc stage, with a collisional cascade having 
been recently initiated (Furlan et al. 2007; Wvatt et al. 2007). The 
orbits of stars are all highly eccentric and inclined, creating a dy- 
namical environment unlike almo st all other known debris discs 
( lTokovinin|[l999l : |Prato et al.l200ll) . 

The dust disc here is generally agreed to be an annulus around 
the B binary, but the exact structure varies from model to model. 
iKoerner et alj ( |2000|) estimate a coplanar narrow ring outwards of 
about 5.0 ± 2.5 au from the two star s, themselves separated by 
approximately 1 au. IPrato et al] ( 1200 ih however determine a 1 au 
high coplanar ring now from 2 and 5 au. lBoden et al.l ( [200 5) argue 
that the line of sight extinction means that th e disc is not copla- 
nar with the sub-binary unless it is very flared. Ipurlan et all ( l2007h 
suggest an inner optically thin ring at 2 au and an outer puffed up 
optically thick wall 0.75 au high at 5.9 au, with a gap between the 
two. Recently, lAkeson et alj ( |2007|) showed that a single continu- 
ous physically and optically thick coplanar disc between 3 and 10 
au can reproduce the observed spectral energy distribution. To ex- 
plain the extinction, they used dynamical models of test particles to 
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show that the inclined stellar orbits could create a warp in the dust 
disc, the outer layers of which could then just intercept the line of 
sight. 

The unusually l arge infrared excess has been argued by 
iLagrange et all ilOOd) to be caused by one of four reasons. The 
first is that a dust avalance is currently in progress. Radiation pres- 
sure acts to push dust grains outwards. When the disc is very dusty, 
these can impact other dust grains, creating more particles that are 
themselves pushed outwards, colliding and creating yet more dust, 
and so on. The second possibility is there is undetected gas main- 
taining the dust population against radiation pressure. The third and 
fourth reasons are applicable when the wide orbit has an eccentric- 
ity of near unity, much higher than currently believed. In this case, 
it is possible that either an encounter with the A pair has heated 
the outer edges of the disc, resulting in abnormally high infrared 
radiation, or disrupted a Kupier Belt-like stucture, resulting in col- 
lisions and releasing large amounts of dust. As the eccentricity is 
now known to be much lower, these two possibilities are less likely, 
although it should be noted that the wide orbit is currently very near 
to periastron. 

The first two possibilities present difficulties in modelling the 
system. If the system is undergoing a dust avalance or is gas-rich, 
then radiation pressur e, collisions and g as dynamics must be in- 
cluded. For example, I Wvatt et al] ( l2007h has calculated the dust 
collision timescale to very short (0.36 years), so collisional effects 
are important in modelling the dust evolution. However, in debris 
discs, dust is generated from an underlying planetesimal popula- 
tion, usually through a collisional cascade. In this case, the dust 
population follows the planetesimal population, which is less com- 
plex to model. If the disc has large amounts of gas or a dust avalance 
is occuring, this may, of course, no longer be the case. Models of 
the planetesimal population remain interesting though, especially 
in such an unusual stellar environment. Indeed, the high inclina- 
tion and eccentricities suggest that little stability is likely, yet some 
must exist. Constraining possible planetesimal locations and then 
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comparing with observations of the dust location may even be able 
to distinguish between the possible scenarios above. They would 
also serve as an indication of whether any planetary stability might 
be possible in this system. 

There is one other debris disc known in a similar stellar system 
to HD 98800. This is GG Tau, which has a circumbinary debris 
disc of several hund red au radius in a qua druple star system with 
the same hierarchy l lGuilloteauet"aLlll999il . However, in this case 
the stars are about ten times more distant, the mass ratios much 
smaller and t heir orbits relativ ely coplanar so more stability would 
be expected teeust & DutrevI [2005 , 2006). Dynamical modelling 
has shown that the circumbinary material forms a sharp ring and a 
more diffuse disc component. 

HD 98800 remains an unusual environment, and the high in- 
clination and eccentricity of the stellar orbits will have a significant 
effect on the dynamics and structure of the planetesimal population 
of the debris disc. This effect has yet to be studied in detail, and 
is investigated here using direct numerical integrations. Section |2] 
gives a description of the stellar system and the simulation param- 
eters. The results are then presented in Section [3] and conclusions 
given in Section]?] 



2 METHOD 

2.1 The Stellar System 

The wide orbit of the sub-binaries A and B is reasonably deter- 
mined, as is that of the stars Ba and Bb in the B stellar pair. How- 
ever the orbit of the other pair, Aa and Ab, is only partly known as 
the smaller star is not resolved. The parameters for all three orbits 
are listed in Table [T] and shown in Figure [T] 

iTokovinin (1999) derives three possible fits to the wide orbit, 
fixing the eccentricity in each case as 0.3, 0.5 and 0.6 and assuming 
a total system mass of 2.6 Mq. Table [T] shows the 0.5 eccentricity 
case, and Table |2] shows the others. Apart from the period (and 
hence separation), the orbital parameters do not vary much between 
the different solutions. 

The A pair is a single-lined spectroscopic binary so only a par - 
tial radial velocity orbit has been determined ^Torres et al.lll995h . 
IPrato et alj j200lh use evolutionary tracks to estimate the mass of 
the Aa star as 1.1 ±0.1 M©. The orbital inclination is not known, 
but the mass and separation of Ab can be found as a function of 
this parameter, as shown in Table [S] A wide range of inclinations 
are still possible even given the constraint that the star is small and 
unobservable. It is likely that this small star is unimportant to the 
dynamics of the circumbinary disc, and so the sub-binary can be 
modelled as a single object. However, it might alter the dynamics 
of the stellar orbits, and this is investigated in the next sub-section. 

2.2 The Stellar Dynamics 

The mutual inclination for the stellar orbits given in Table[T]can be 
calculated using the formula 

cos (j) = cos ii cos 42 + sin ii sin 12 cos(r22 — ^i) (1) 

where fij are the longitudes of ascending node, ij are the inclina- 
tions relative to a given reference plane and (j) is the mutual inclina- 
tion, the angle between the angular momentum vectors of the two 
orbits (see e.g. ISmar^ll953h . It is found to be in this case 143.7°, 
and is similar for the other two orbits (see Table [2l(. Notably, it is 
retrograde. 



Table 1. The orbital and physical para meters for t he ste llar system HD 
98800. The wide orbit of AB is fit II from 'Tokoviniiil il999l) . with the semi- 
major axis calculated assuming a total system mass of 2.6 Mq and the MID 
taken from the middle of the year 2025 given as the time of pe riastron of the 
AB pair The Bab pair orbit is taken from the joint-fit given bv|Boden et alj 
1 2005). The orbit of the A sub-binary is from lTorres et al.l i 19951) . Note that 
the reference plane is that perpendicular to the line of sight and the longi- 
tudes of the ascending nodes ai'e measured from North thi'ough East. 



Orbital Wide A sub-binary B sub-binaiy 

Parameter A Aa Ab Ba Bb 



Mass(MQ) - M±0.1 0.699±0.064 0.582±0.051 

Radius (Rq) - 1.09 ± 0.14 0.85 ±0.1 1 

a (au) 67.6 0.447±0.0I3 0.536±0.013 

Period 345 yr 262.15±0.51 d 314.327±0.028 d 

e 0.5 0484±0.020 0.784t)±0.0053 

88.3 66.8±32 

uiC) 224.6 644±2.1 289.6±l.l 

nC) 184.8 337.6±24 

t(MJD) 60840 8737.1±1.6 52481.34±0.028 



Table 2. The three different orbital fits for the wide orbit AB from 
TokovinirJ il999l) . Here <j) is the mutual inclination to the orbit of the B 
stellar pair 



Orbital 
Parameter 


I 


Orbit 
II 


III 


a (au) 


61.9 


67.6 


78.6 


Period (years) 


302 


345 


429 


e 


0.3 


0.5 


0.6 




87.4 


88.3 


88.7 


"(°) 


210.7 


224.6 


224.0 


no 


184.8 


184.8 


184.8 


T (MJD) 


59379 


60840 


61205 


4>n 


143.0 


143.7 


143.9 



As mentioned, it is likely for an investigation of the dynamics 
of the circumbinary disc that the Aab system can be reasonably 
approximated by a single star of mass 1.3 Mq, consistent with 
the minimum ma ss of Ab and the total system mass assumed by 
iTokoviniiJ ( Il999h . In this case, the system becomes a hierarchi- 
cal triple and can be n umerically studied using the MoiRAI code 
l l Verrier & Evansl2007t) . A million year simulation of orbit II in this 
case is shown in Figure |2] Energy is conserved to a relative error 
of 10~^ and the integration was ch ecked with a standard Bulirsch- 
Stoer integrator jPress et ai]|l989h and found to be in agreement. 
The semimajor axes of the stars and the wide orbit's eccentricity 
are not shown, as they remain constant through out the simulation. 
Interestingly, the system is currently near its maximum eccentricity 
and mutual inclination. For all three possible outer orbits, the ec- 
centricity and mutual inclination vary smoothly in a secular man- 
ner, with maximum angular separation between the orbital planes 
occurring for minimum eccentricity of the inner orbit. The ampli- 
tude of the variation is very similar in all cases and the periods are 
65, 60 and 75 Kyr for orbits I, II and III respectively. In fact the sys- 
tem remains in this stable configuration for at least a Gyr, as well 
as to at least 10 Myr ago. 

The orbital beh aviour is well described by the octupole secular 
theory o f lFordetalJ ( l200Q) , as overplotted in Figure|2] This theory 
uses a third order expansion of the systems Hamiltonian in the ratio 

Table 3. The mass of star Ab and the A binary's semimajor axes as a 
function of inclination to the plane of the sky. 

Inclination (°) 90 70 60 50 40 30 20 10 

MAb (Mq) 0.22 0.23 0.25 0.29 0.36 0.49 0.80 2.36 
a(au) 0.88 0.88 0.89 0.89 0.91 0.94 0.99 1.21 
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Figure 1. The physical orbit of the quadruple star system HD 98800 in 
three dimensions. The wide orbit is fit II from Tokovinin ( 1999) and the B 
pair orbit that of Boden et al., ( 2005). Units are in au, and the plot is centred 
on the barycentre of the B stellar pair. A projection of the orbits onto the 
plane of the sky is also plotted along with the relevant directions. The B 
parr orbit is overplotted magnified 20 times as dotted fines. The position of 
stars at the epoch MJD 52481.34 (mid 2002, inner orbit at periastron) and 
the direction of the orbits are given. Finally, the line of intersection of the 
two orbits is also plotted as a dashed line. 



of the semimajor axes of the stellar orbits to obtain a set of coupled 
first order differential equations, their equations (29) to (32), for 
the time variation of the eccentricities and arguments of periastron 
of the inner and outer orbits that can be numerically solved. The 
inchnation and nodes of the system are then derivable from these 
elements, and the semimajor axes remain constant. From the fig- 
ure, it can be seen that the theory is in excellent agreement with the 
results from the full equations of motion, with only a small phase 
drift after a Myr. In fact, because the two orbits are fairly sepa- 
rated a quadrupole level theory is sufficient to describe the system. 
In this approximation the outer eccentricity is a constant, as seen 
here. Thus, in the three body approximation the stars follow a stable 
secular evolution and can be accurately integrated by the MoiRAI 
code. 

The full four body system can now be considered. The mass 
of Ab and the semimajor axis of the orbit can only be resolved 
by assuming the inclination to the line of sight (see Table [3]l. The 
longitude of ascending node of the orbit remains an unknown, but 
is needed to constrain the mutual inclination of the orbital planes. 
Therefore, to investigate possible dynamics, a set of simulations 
were run for inclinations in the range 90° to 30° (shown in Table[3} 
with values of ascending node ranging from 0° to 315° in 45° steps 
and using wide orbit II. To do this, the fourth star was approximated 
as a planet around its primary, Aa, in the MoiRAI code. 

In all cases, there is little difference to the three body results. 
An example is shown in Figure [3] The only change is to the secu- 
lar period of orbit Ba-Bb together with a slight modulation of their 
minimum eccentricity, and in some cases even this does not oc- 
cur. A slight change in the period of the secular variations is un- 
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Figure 2. The three body stellar orbital evolution for orbit II over 1 Myr 
from numerical and theoretical modelling. The blue line sfiows the results 
from the symplectic integrator, identical to those from the Bulirscfi-Stoer 
shown in red. The octopole theory results are shown in black, and differ 
only slightly in phase. The semimajor axes of both orbits and eccentricity 
of the outer orbit are constant and not shown. 
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Figure 3. The evolution of the stellar orbits for the four body case with 
*Ai) = 30° (relative to the plane of the sky) and QAb = 45° compared 
to the three body approximation. The left hand panels show the orbit of the 
inner binary Ba-Bb, red indicating the three body approximation and black 
the four body case. The right panels show the evolution of the orbit of the 
Aa-Ab binary. The eccentricity of the wide binary A-B once again remains 
effectively constant, as do the semimajor axes of all orbits. 



likely to alter the overall st ructure of the planetesimal disc. In- 
deed, iBeust & DutrevI ( l2006h also find that modelling the distant 
sub-binary in GG Tau as a single object has little effect on the disc 
structure there. Hence, the three body approximation is a reason- 
able assumption and will be used for the purposes of this investiga- 
tion. 



2.3 Other Physics 

As the disc mass is low (0.002Mm. |Prato et al.ll200ll) . It is reason- 
able to model the planetesimals as non-interacting test particles. 
However, there is one interaction process that may be important in 
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the modelling of the planetesimals, namely collisions. Following 
IWvatt et al] ( l2007h . it is possible to esimate the collision timescale 
from the observed infrared excess. Although this assumes a colli- 
sional cascade is underway, it should provide a lower limit to the 
timescale if dust is in fact maintained by gas or generated in an 
avalance. The collision timescale is given by their equation (13) as 



Table 4. The test particle grid u.sed in the simulations. 



tc = 



3.8pr^-^drDc 
M°-^Mtot 



/(12g_ 



20)(1 + 1.25(e/7)^ 



(18-9g)G(g,X,) 



(2) 



where q is assumed to be 11/6 for a coUisional cascade and G is a 
function of q and Xc ~ Dcc/Dc, Dc being the diameter of largest 
planetesimal in the cascade, taken as 2000 km, and Dec being the 
smallest planetesimal that has enough energy to destroy another of 
size Dc- p is the planetesimal density and taken as 2700 kg m^^, 
r is the planetesimal belt radius in au, dr is the planetesimal belt 
width in au, Mtot is the total mass of material in the cascade in 
Mfg, e is the mean orbital eccentricity of planetesimals and I is 
the mean orbital inclination of planetesimals in radians. The factor 
G{q, Xc) is given by their equation (9) as 

6g- 10, 



G{q,Xc) 



3g 



3q 



+ 



-1)- 



3(7-3 

and Xc by their equation (1 1) as 



3q~^ 
39 _ 1) 



-ix. 



4-3ij 



1) 



Xc 



1.3x10-^ (Qj,rACV(e,/)"') 



2\l/3 



(3) 



(4) 



where /(e, I) is the ratio of the relative collision velocity to the Ke- 
plerian velocity and taken as 0.1 and is the dispersal threshold 
which is the specific incident energy needed to destroy a particle, 
assumed to be 200 J kg^ ^ . In addition, it is also assumed that e ~ /. 

The two stars Ba and Bb are approximated as a single ob- 
ject of mass and l uminosity as 1.57 Mq and 0.58 Lq respectively 
jPrato et al.l200ll) . The total mass of the disc can be calculated from 
equations (4) to (6) of Wyatt et aL t2007.) as 



Mtot = 5.194o-tot = 14.36r 

and hence the collision timescale as a function of r and dr is 

1 



tc = 1.014 X 10V°-^ 



dr- 



(5) 



(6) 



Giq,Xc) 

The exact extent of the planetesimal disc is as yet unknown, but the 
location of the dust can be used as an approximation. Th e minimum 
radius for dust to be able to exist is given as 2.2 au by Low et al. 
( I2OO5I) , but lKoerner et al. ( 2000) estimates the disc to be between 
5.0 and 18 au. In the first case of a wide disc around 2 au, the col- 
lision timescale is 30 Kyr, and in the second case of a disc from 
5 to 18 au, it is 200 Kyr. Preliminary results from test simulations 
indicated that the disc must lie further out than 2 au, so the higher 
timescale is applicable and gravitational effects will dominate the 
behaviour of the planetesimals here. Thus, for the purpose of de- 
termining overall disc structure, no additional physics needs to be 
included. 



2.4 Method of investigation 

To summarise, the stars are modelled as a three body system with 
a disc of massless planetesimals interacting through gravitational 
effects with the stars only. The three possible orbital configurations 
for the wide binary will all be considered. To model the disc, the 
MoiRAi code has been shown to be accurate, and the test particle 
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Figure 4. The test pailicle decay rates for the thi'ee simulations. The half- 
life (ignoring the first 0.1 Myr when all initially unstable particles are re- 
moved) is about 0.3 Myr for orbit I, 0.2 Myr for orbit II and 0.1 Myr for 
orbit III. 



disc can be implemented as circumbinary particles around the B 
stellar pair. 

The models discussed in the introduction have all assumed 
that the disc is coplanar with the B binary, but there is evidence 
for inclined discs around similar multiple stars. For example, the 
close binary pair in the T Tauri system were recently determ ined 
to have misaligned circumstellar discs jSkemereta"l]|2007h . and 
polarization surve ys have found a small number of similar cases 
jMonin et alj2006 ). Because of this, the test particle distribution is 
not initially taken as coplanar. Instead, they are spaced uniformly 
in inclination and longitude relative to the B binary pair. The initial 
surface density of the disc is taken as 1/r and the grid of test par- 
ticles runs from 1.75 au to 33.85 au (inner binaries radius to half 
the outer binaries radius), as shown in Table|4] As the present day 
configuration, stability and geometry of the system are of interest to 
compare to observations, the simulations are run from 1 Myr ago to 
present day (defined as the periastron time of the Ba-Bb pair from 
Table [T). The results from these simulations are presented in the 
next section. 



3 RESULTS 

The stability of planetesimals is indicated by the stability of the test 
particles. These are removed from the simulation if they cross either 
of the stellar orbits or if they become unbound. The test particle de- 
cay rates are shown in Figure|4]for the three different simulations. 
The half-life is short and by the end of the simulations the decay 
rate has levelled out, although some particles are still being slowly 
eroded from very unstable locations. The majority of unstable par- 
ticles are removed on 10 to 100 Kyr timescales as they become 
unbound or cross the outer orbit. If the simulations are run for an 
additional million years, there is no further significant loss of par- 
ticles so the simulation length is sufficient to determine the system 
state, especially given the system age of about 10 Myr. 

The final structure at the end of the simulation is similar in 
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Figure 5. Edge and face on views of the circumbinai'y material for tlie simulation using orbit II. The inner binary's orbit lies in the 2 = plane and the line of 
intersection with the outer stellar orbit is along the y-axis. The orbits of the stars and their current positions are overplotted in red, but the wide orbit is shown 
at a tenth of its actual size. The line of sight is also shown as a dashed line. 
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Figure 6. The final test particle distributions for the simulation using orbit II. The two panels on the left show the incHnation relative to the inner binary as a 
function of semimajor axis and radial distance from the center of the binary, with colour indicating inner mutual inclination as shown. The two panels on the 
right show the inclination relative to the wide orbit, but the colour still indicates inner inclination. Oveiplotted as a dashed line is the mutual inclination of the 
two stellar orbits. 



all three orbital cases, and is well illustrated by the case of orbit 
11. Shown in Figure|5]is the spatial distribution of the circumbinary 
material. A warped coplanar disc is apparent around the inner bi- 
nary, as well as a large amount of higher inclination material. The 
disc appears to have a small inner ring separate to the main bulk 
of particles. There is also a slight clumping of material apparent, 
perpendicular to the inner orbits line of apses. The high inclination 
material appears to form another ring or halo perpendicular to the 
disc, most obvious in the middle plot of the figure. 

These different structures are more clearly seen in a plot of 



inclination as a function of distance from the inner binary, as shown 
in Figure |6] Here, the particle distribution relative to the iimer and 
outer orbits is plotted for both radial distance and semimajor axis. 
To identify different populations, particles are colour coded to their 
inclination relative to the inner binary's orbit. 

Three distinct features are seen in these plots: a prograde 
coplanar disc with two gaps; a retrograde disc with one gap and 
an extended ring or halo, as seen in the spatial plot, from 50° to 
130° inner inclination. These populations are well separated in in- 
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clination by two gaps around 50° and 140°, the later being centred 
on the inchnation of the outer binary's orbit. 

The coplanar disc extends from about 4 to 15 au, with two 
gaps between 4.5 and 5 au (seen very clearly in Fig.[5]l and 7 and 8 
au. The gaps are less obvious in the radial distribution due to parti- 
cle eccentricity, but are still apparent. The outer regions of this disc 
are inclined up to almost 50° and provide the warped material seen 
in the spatial plot. There is also a series of breaks in the semimajor 
axis distribution slightly apparent both here and in the retrograde 
disc (these are very clear in Figure[9j, presumably due to resonant 
features as they are much larger than the grid resolution of 0. 1 au. 
Relative to the outer orbit, the disc material in the outer warped 
regions is perturbed up to almost completely retrograde orbits, ex- 
plaining its stability. 

The retrograde (relative to the inner binary) disc is smaller, 
extending from 4 to 9 au, but also has a gap between 4.5 and 5 
au. The outer regions here are perturbed towards the same inclina- 
tion as the outer orbit, similar but opposite to the coplanar disc's 
structure. This disc is in fact much sharper relative to the outer star, 
where it is prograde with an inclination of around 30° . 

The last component is the halo-like structure. This extends 
only from 5.0 to 8.0 au but covers a large range of inclinations 
relative to both stellar orbits. This material is unusual in remaining 
stable at very high inclinations, b ut similar stable high in clination 
particles have been seen before bv lBeust & DutrevI bOOeh in simu- 
lations of GG Tau. 

These three separate populations are also very distinct in a plot 
of final versus initial inclination, as shown in Figure |7] Particles in 
the prograde disc start with initial inner inclinations in the range 
0° to 50° and remain within that range. Those in the halo, starting 
in the range 60° to 125°, also remain there although there is less 
variation in the middle of this group. Finally, the retrograde disc 
is clearly confined by the inclination of the outer stellar orbit, and 
particles here start and remain in the range 145° to 180°. These 
populations are again clear in the outer inclination distribution, as 
can be seen from the colour coding of the particles. There is in 
fact a fourth diffuse population visible here, at around 70° to 90°, 
which forms a more diffuse ring at a different angle to that of the 
main halo. 

The initial outer inclination distribution is different for the or- 
bit I and III simulations, as the stars start in slightly different places. 
However, the components evolve to the same final inclination dis- 
tribution, showing that the results are robust and do not strongly 
depend on parameters such as the initial longitudes. Note also that 
in the orbit II case no particles start with initial outer inclinations 
less than about 30°, so it could be possible that another region of 
stability exists here. However, a disc of particles started coplanar 
with this outer orbit very rapidly became unstable, with no material 
surviving to the end of the simulation. 

A radial profile of the three main structures can be plotted 
by using the inner inclination range to characterise them. This is 
shown in Figure [8] for each simulation. The coplanar disc (black) 
is defined as starting between 0° to 50°, the halo (green) between 
55° to 140° and the retrograde disc between 145° to 180°. In each 
case the profiles of each component are very different, and the gaps 
in the two disc very apparent. 

The radial profiles here are the result of the evolution of an ini- 
tial 1/r distribution. A simulation was rerun for a flat initial density 
profile in the orbit II case. The resulting disc profiles were similar, 
although the inner and outer edges were slightly steeper. The profile 
of the halo also remained largely unchanged. 

The most notable difference between the three simulations is 




Initial inner mutuol Inclination {") Initial outer mutual Inclination {") 



Figure 7. The initial test particle inclinations related to their iinal incli- 
nations for the simulation using orbit II. The left panel shows inclination 
relative to the inner binary and the right panel relative to the outer wide 
orbit. Colour now indicates initial outer inchnation to highlight the small 
population around 80° in the right hand plot. The dotted vertical line shows 
the initial mutual stellar inclination and the horizontal the final. 

the generally greater stability in the lower eccentricity cases and 
the existence of an extended retrograde disc in the orbit I case. This 
stability is shown in greater detail in Figure |9l and compared to 
the other simulations. As the extended prograde coplanar disc is 
retrograde relative to the outer star, this disc is prograde to it. The 
extended prograde disc is itself less populated here, but both fea- 
tures can partly be explained by the initial particle inclinations. The 
stellar orbits are different here and particles that are retrograde rel- 
ative to the inner star are more coplanar with the outer star than in 
the other orbital cases, and those that have prograde inner inclina- 
tions have lower outer inclinations. Thus, fewer particles start in the 
prograde extended disc, while more start in the retrograde case. It 
should be noted, however, that even if the initial inclination distri- 
bution is similar to the other two simulations, a retrograde extended 
disc is still seen. The eccentricity of the stellar orbits is therefore 
still likely to be important to these features, as can be seen by the 
decreasing radial extent of the prograde disc from orbit I to III in 
Figure |9] Indeed, as discussed above, a retrograde disc relative to 
the outer star for the orbit II case is completely unstable. 

The radial extent of the structures seen in the three simula- 
tions is detailed in Table |5] and compared to both previous esti- 
mates of the size of the observed dust disc and empirical stability 
limits from numerical studies of general binary and triple systems. 
These empirical limits show the inner and outer radius for coplanar 
circumbinary s tability by mode lling the stars as a low inclination 
triple system d Verrier & Evans 2007) and as two coplanar decou- 
pled bina ry systems ^Holman & Wiegert. 1999.) . A more general re- 
sult from lVerrier & Evanj ( l2007h is that in cases of high stellar ec- 
centricity and mass ratio the stable region was more complex than 
a completely stable ring, with gaps appearing that were most likely 
due to overlapping resonances from the two stellar orbits, as seen 
here. The extended outer regions of the prograde and retrograde 
discs are outside the empirically predicted stable region. This is 
most likely an effect due to the high inclinations and the particles 
retrograde nature relative to the outer or inner orbit, as retrograde 
orbits are generally more stable than prograde. However, the inner 
edges are well predicted. 

An important feature of the observations of the system is the 
extinction and photometric variability towards the B binary, at- 
trib uted to disc material along the line of sigh t by Tokovinin (199^ 
and lBoden et"ai] ( l2005l) . lAkeson et alj ( l2007h find that for the sys- 
tem they investigate the line of sight can just intercept the top of a 
warped disc of test particles. To see what material, if any, occurs 
along the line of sight here, we plot the azimuthal structure of the 
disc, as shown in Figures [TO] and [TT] for the case of orbit II. Here, 
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Figure 8. The radial profile of the three inclination components of the circumbinary material for all three simulations, showing the surface density as a 
function of radius. The coplanar disc is plotted in black, the halo in green and the retrograde disc in red. Note that for the orbit I case the retrograde disc starts 
at 105° instead of 145° due to the differences in the initial inclinations of the stars for this orbit. The retrograde disc and halo here also now overlap in inner 
inclination, and are separated using their outer inclination instead. 



Table 5. The radial extent of the disc around HD 98800 B. The top four 
rows show other models of this system, the rest show empirical fits from 
general binary and triple systems and the extent of the thi'ee different incli- 
nation components for each possible outer orbit from the simulations here. 



Model Inner edge (an) Onter edge (an) 



Koerner et al. (2000) 
Praloelal. (2001) 
Fnrlan et al. (2007) 
Akeson et al. (2007) 


2.0 
2.0 
3.0 




5.0 ± 2.5 
gap 








5.0 
5.9 
10.0 


Oibil I 


Holman & Wieaert (1999) 
Verrier & Evans (2007) 


4.1 

3.9 












10.9 
11.0 


Disc 


3.0 


-4.0 


gap 4.5 - 


8.0 


gap 


9.0- 


20.0 


Halo 


5.0 












8.0 


Retrograde 


3.5 


-4.0 


gap 4.5- 


7.0 


gap 


7.5 - 


13.0 


Orbit 11 


Holman & Wie.gert (1999) 


4.1 












4.3 


Verrier & Evans (2007) 


3.9 












0.6 


Disc 


4.0 


-4.5 


gap 5.0- 


7.0 


gap 


8.0- 


15.0 


Halo 


5.0 












8.0 


Retrograde 


4.0 


-4.5 


gap 5.0 








9.0 


Orbit III 


Holman & Wieeert (1999) 
Verrier & Evans (2007) 


4.1 

3.9 












7.0 
7.3 


Di.sc 


3.0 


-3.5 gap 4.0- 


-4.5 gap 5.0- 


7.0 


gap 


8.0- 


15.0 


Halo 


5.0 












8.0 


Retrograde 


4.0 


-4.5 


gap 5.0 








8.0 



the system has been rotated so that the inner binary lies in the hor- 
izontal plane with its periastron at an azimuthal angle of 0° . Each 
panel then shows the height above the plane as a function of radial 
distance within the plane of the inner binary in segments of 10° in 
azimuthal angle. Particles are colour coded to inner mutual incli- 
nation as for Figure|6]to distinguish the different components. The 
locations of the three stars are also shown (the position of the outer 
star is shown reduced by a factor of ten), as is the line of sight on 
the relevant plot. 

The warp in the prograde disc is very apparent, with material 
here reaching heights of almost 10 au above the plane. The retro- 
grade disc is slightly less warped, but does not extend out as far. The 
halo material shows up very clearly as a ring, instead of a contin- 
uous shell covering all angles. There is a small amount of material 



perpendicular to the main ring only at very high z, the fourth pop- 
ulation visible in Figure|7] The warp lies along the 30°-210° line, 
as illustrated in the lefthand panel in Figure[TT] but the line of sight 
at around 160° still intercepts a large amount of perturbed copla- 
nar disc material, a good 5 au above the plane of the binary's orbit. 
Figure[T2]shows similar plots for the other two orbital cases for the 
azimuthal segment containing the line of sight. In the high eccen- 
tricity case of orbit III, almost no material is intercepted, while in 
the low eccentricity orbit I simulation far more prograde and retro- 
grade material remains, making either this or orbit II the most likely 
orbital configurations if a warped planetesimal disc and associated 
dust is the cause of the observed extinction. 

Over the length of the simulation, the warp processes with a 
timescale equal to twice that of the secular period of the stars, fol- 
lowing almost perpendicular to the circulation of the line of inter- 
section of the two stellar orbits. The extent of the warp in fact de- 
creases as the mutual inclination between the stars increases (since 
the orbit is retrograde the higher the inclination the closer the two 
planes). The non-symmetrical distribution seen in Figure [TT] per- 
sists, with one side usually greater in height than the other. The 
warp is not a short term feature, and the system will persist in its 
current configuration for some time, so if dust follows the planetes- 
imal distribution it is very likely the extinction is indeed caused by 
the warp in the disc. It should also be noted that the warped material 
remains after an additional Myr, and is not slowly eroded away. 

The azimuthal particle distribution shown in the right hand 
panel of Figure [TT] shows that the two discs are not skewed in any 
particular direction, but that the halo is aligned along the minor 
axis of the inner binary's orbit. In fact, it remains aligned at this 
angle over the entire simulation. The material in this halo is seen 
in projection in the right most panel of Figure|5]as the two clumps 
discussed earlier. As this material follows the pericentre of the inner 
orbit, which processes on the secular time scale, it will not intercept 
the line of sight for some time and is unlikely to be related to the 
observed variability and extinction. 

An important final question raised bv| Pratoetal] ( l200lh is the 
lack of a similar circumbinary disc around the other binary A. As 
the eccentricity and mass ratio of this pair are both smaller, a disc 
should be more likely here. In fact, the empirical criteria place the 
stable zone from around 3 au to II, 8 and 7 au for orbits I, II and III 
respectively. A preliminary simulation taking A as a single star con- 
firms this greater stability, so the question still remains as to why 
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Figure 11. The azimuthal distribution of particles for tlie simulation using 
orbit 11. The left hand panel shows a polar plot of the waip in the prograde 
(black) and retrograde (red) discs. The halo is not a disc with respect to this 
plane so not plotted. The average height away from the plane of the inner 
binary is shown as a function of azimuthal angle, increasing anticlockwise 
from the horizontal axis, and the periastron of the inner binary is at 0° . The 
line of sight is shown as a solid line and the line of intersection of the two 
orbits as a dashed line (the rising node is towards the bottom right). A plus 
marks the side of the warp that is above the plane. The right hand panel 
shows a polar plot of the number of particles in each 10° angle bin. Red 
and black are as before and the halo is now also plotted in green. 



there is no observed disc here. It is possible that the inclination of 
this binary's orbit places it so that no planetesimals can remain sta- 
ble, and if so this would provide limits on the orbit of this stellar 
pair. Another alternative if there is indeed a disc of stable planetes- 
imals here is that it is far less dusty and so unobserved. 



4 CONCLUSION 

Dynamical simulations of a planetesimal population in the debris 
disc around the Bab stellar pair in HD 98800 have been run. By 
studying a wide range of inclinations, three distinct stable popula- 
tions have been identified. These are a prograde disc, a retrograde 
disc and a high inclination halo. The radial profiles of each com- 
ponent are different and distinct. The discs both have large radial 
gaps, caused presumably by overlapping resonances from the stel- 
lar orbits. The radial extent of the discs are summarised in Table |5] 
but are generally from 3 to 15 au for all three orbits, with gaps at 
around 4.5 au and 8 au. 

The line of sight can currently pass through slightly warped 
material in the prograde disc, and this would account for the ob- 
served extinction and variability if the dust distribution followed 
the planetesimals. However, this alignment with the line of sight 
effectively only occurs for the I and II orbits (with eccentricities of 
0.3 and 0.5 respectively), which would rule out the higher eccen- 
tricity case (e = 0.6) as a possible orbit for the outer star 

The radial profile is complex, and illustrates that, for discs in 
multiple systems, dynamical effects purely due to the stars are very 
important. Indeed, if most triple stars are within the limit for reso- 
nance overlap to sculpt the stability zone, then this has large conse- 
quences for debris discs and planet formation in such systems. 

The simulation results here can be compared to other mod- 
els of the system and estimates of the dust distribution. Indeed, 
as previously mentioned, comparisons to observations of the dust 
could place constraints on the physical processes occuring in the 
circumbinary disc. The bulk of the dynamically stable planetesi- 
mals are between 5 and 7 au. This matches up well with the pre- 
diction of Koemer et alj i2000t) of a ring outwards from 5 au. The 
model of IPrato et alj ^2001 ^) places the disc from 2 to 5 au, which, 
apart from the ring around 4.5 au, is unstable here. iFurlan et al.l 



feOOTl) suggest an inner ring at 2 au and a thicker puffed up com- 
ponent at 5.9 au, the only model to predict gaps and similar to the 
prograde discs here. The later of these components matches up to 
the bulk of the material here, and there is in deed an inner ring seen, 
just further out at 4 au. IPrato et alj (1200 ih estimate the height of 
the disc as 1 au, and lFurlan et al.l ( l2007h as 0.75 au. However, here 
material in the prograde planetesimal disc is not only warped but 
very flared, reaching heights of 5 to 10 au. This is an important dy- 
namical feature that should be taken into account in models of the 
observations. Although roughly agreeing in location, these models 
do not have enough detail as yet to further compare to the planetes- 
imals. 

The dynamical model of lAkeson et alj ( l2007h finds stability 
from 3 to 10 au, which does not match up that well to the limits 
here, and in fact their resulting disc has a very different geometry. 
Their model uses different orbital parameters for the stellar system, 
most notably a different mutual inclination - as their aim is not 
to reproduce the current configuration but to look at the general 
morphological structure of debris discs in highly inclined triple star 
systems. They also find that the line of sight can intercept a warp in 
the disc, but due to the different geometry it crosses at a different 
point and at the maximum warp in the disc, whereas here the warp 
is only marginally orientated towards the line of sight. It is in fact 
material in the sparser populated outer regions of the prograde disc 
that intercepts the line of sight here, from 8 to 15 au, which is a 
region that in a coplanar stellar system is predicted to be unstable. 

The possible dynamical reason for the high infrared excess al- 
ready mentioned in the Introduction is a close pericentre passage 
of the A binary stiring up planetesimals and resulting in high colli- 
sion rates. Since the binary's period is about 300 years, many such 
pericentre passages have occured, and by the end of the simulations 
the particles appear very stable with regards to this. For example, 
there appears to be no periodic increase in planetesimal eccentric- 
ities as the outer star passes close to the disc. There is, however, 
another feature that may cause higher than normal collision rates. 
The extent of the flare and warp in the extended disc decreases as 
the acute angle between the stellar orbits decreases (and the mu- 
tual inclination increases). As this occurs, planetesimals that have 
been spread out in inclination are now packed into less space. This 
is likely to result in an increased number of collisions and is par- 
ticularly relevant as the current stellar configuration is such that 
mutual inclination is large and the disc is near its narrowest point. 
However, modelling including collisions is needed to quantify this 
effect. 

The stable high inclination particles in the hal o are curious 
from a dynamical perspective. The Kozai mechanism ( lKozaill963) 
might be expected to remove such particles very quickly - however 
these particles do not seem subject to this. It is possible that they 
are in fact an artifact of the numer ical method, but checks with a 
standard Bulirsch-Stoer integrator jPress et alj|l989l) and the sim- 
ilar observation in simulations of GG Tau discount this. The high 
stellar inclinations involved in the system may be one explanation, 
but it is worth looking at the stability should one of the binaries be 
removed. If the outer star is not present, a sharp inner disc edge is 
seen at around 4 au for all inclinations, with all particles outwards 
from this remaining stable. If the inner binary is approximated as 
a single star instead, then a reasonably sharp outer edge is seen at 
around 7 to 10 au, for the lower inclinations only. The near polar 
inclinations in this case are now unstable. The high inclination par- 
ticles are within this region which must be shaped by the combined 
perturbations from both stellar orbits, so are in a new regime of 
dynamical behaviour. The mechanism stabilising these particles, as 
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Figure 9. Stability maps comparing the results from the three simulations using orbits 1, 11 and 111. The average survival time is shown as a function of initial 
semimajor axis and inner inclination, with black indicating all particles starting at a given grid point remained at the end of the simulation through to white 
indicating the location was very quickly unstable. 



well as that stabilising the extended prograde and retrograde discs, 
is the subject of a future paper. 

Some consideration is needed, however, of how such material 
would form. T Tauri multiple systems are common and believed to 
form primordially through fragmentatio n processes, wh ich could 
result i n non-coplanar disc s ( Bate.2000 , Bate et al.l200oh . As men- 
tioned, |Mo^net^] j2006l) find a small number of multiple T Tauri 
systems with misaligned discs in their polarization survey, but sug- 
gest that perhaps these are perturbed disc soon to be realigned. So 
there is some evidence that non-coplanar particle distributions are 
plausible. It may be possible as well that particles from the disc can 
be captured into the polar stability region, and certainly a coplanar 
particle distribution is quickly perturbed to fairly high inclinations. 

There is much further work to be done in modelling this sys- 
tem. A model of the planetesimals including collisions, although 
expected to be a minor effect on the dynamics, would indicate if the 
rate was currently high and may explain the unually dusty nature 
of the system. A model that also included dust and dust collisions, 
as well as interactions with any gas in the system, would then fully 
model the system and allow detailed comparisons with the obser- 
vations. Although not an easy task given the nature of the system 
it would narrow down the physics and dynamics at work, impor- 
tant to our understanding of this stage of the planetary formation 
process in multiple star systems. 
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Figure 10. The azimuthal distribution of material for tlie simulation using orbit 11. The panels show the radial distance within the plane of the inner binary 
and the height above it in segments of 10°. The periastron of the inner binary is at 0°. The stellar positions are overplotted with the outer stai" shown reduced 
by a factor of ten, and the line of sight shown as a dashed line. The test pailicle colour indicates current inner inclination, as in Figure[6l 
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Figure 12. The material intercepting the line of sight in all three simulations. As for Figure ffol radial distance within the plane and height above it is shown 
with plot points colour-coded to initial inclination. Here however only the segment near the line of sight has been shown for each simulation. 
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